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ABSTRACT 

We combine a high-resolution hydro-simulation of the ACDM cosmology with two radiative transfer 
schemes (for continuum and line radiation) to predict the properties, spectra and spatial distribution of 
fluorescent Lya emission at z ~ 3. We focus on line radiation produced by recombinations in the dense 
intergalactic medium ionized by UV photons. In particular, we consider both a uniform background 
and the case where gas clouds are illuminated by a nearby quasar. We find that the emission from 
optically thick regions is substantially less than predicted from the widely used static, plane-parallel 
model. The effects induced by a realistic velocity field and by the complex geometric structure of the 
emitting regions are discussed in detail. We make predictions for the expected brightness and size 
distributions of the fluorescent sources. Our results account for recent null detections and can be used 
to plan new observational campaigns both in the field (to measure the intensity of the diffuse UV 
background) and in the proximity of bright quasars (to understand the origin of high colum-density 
absorbers). 

Subject headings: cosmology: theory - intergalactic medium - large-scale structure of universe - line: 
formation - quasars: absorption lines - radiative transfer 



1. INTRODUCTION 

Hydrogen absorption-line systems observed shortward 
of Lya emission in quasar spectra constitute an im- 
portant probe of the physical state of the intergalac- 
tic medium at high-redshift. These spectral features 
are shaped by the combined action of gravity, hydrody- 
namics and photoionization processes which determine 
the local density and the velocity field of neutral hydro- 
gen within the absorbers. Numerical simulations suggest 
that the so called Lyman-a forest is generated by dif- 
fuse, sheetlike and filamentary structures with a mean 
density which is between 1 and 10 times higher than 
the cosmic average (Cen et al. 1994; Zhang, Anninos 
& Norman 1995; Hernquist et al. 1996; Miralda-Escude 
et al. 1996). These low-column-density systems are 
highly ionized by the extragalactic background of Ly- 
man continuum photons generated by young stellar pop- 
ulations and quasars. At the opposite extreme, Lyman- 
limit (LLS, A'hi > 10 172 cm -2 ) and damped Lyman-a 
(DLA, TVhi > 10 20 ' 3 cm~ 2 ) systems correspond to con- 
centrations of atomic hydrogen which are optically thick 
to the cosmic ionizing background. Numerical simula- 
tions suggest that they arise in dense gas clouds with a 
meatball topology. On cosmological scales, they appear 
to form a collection of isolated clouds which trace the 
cosmic web. 

Optically thick clouds are expected to emit fluores- 
cent Lya photons produced in hydrogen recombinations 
(Hogan & Weymann 1987; Gould & Weinberg 1996). 
This emission is concentrated in the outer parts of the 
clouds where hydrogen is significantly ionized by the ex- 
ternal UV background (tll ~ 1). However, Lya photons 
cannot directly escape the clouds because of the large 
optical depth in the center of the line (TL yQ ~ 10 4 tll)- 
Each photon thus suffers a large number of resonant scat- 
terings (more precisely: absorptions and re-emissions) by 
neutral hydrogen atoms in the ground state. Each scat- 



tering adds a small Doppler shift to the frequencies of 
the photons due to the thermal (and turbulent) motions 
of the atoms. Therefore, photons execute a random walk 
both in frequency and in physical space until their fre- 
quencies are shifted sufficiently away from the line center 
and they are able to escape the medium in a single flight 
(Zanstra 1949). 

Monte Carlo simulation (e.g. Ahn, Lee & Lee 2001; 
Zheng & Miralda-Escude 2002b and references therein) 
is the most popular method for addressing the radia- 
tive transfer problem. Analytical solutions only exist for 
highly symmetric systems. For instance, the emerging 
spectrum from a plane-parallel and static homogeneous 
slab is characterized by two sharp peaks in the Doppler 
wings of the line (Neufeld 1990 and references therein). 
The plane-parallel solution approximately holds also for 
self-shielded systems where the ionized layer which sur- 
rounds the neutral region is thin with respect to the char- 
acteristic radius of the cloud. In this ideal case, opti- 
cally thick systems act as efficient mirrors which convert 
nearly 60% of the impinging ionizing flux into Lya pho- 
tons (Gould & Weinberg 1996). 

Direct imaging of fluorescent sources would lead to a 
major advance in our understanding of galaxy formation. 
Determining the size distribution of LLS at z > 3 would 
be crucial to distinguish whether they arise from pho- 
toionized clouds in galactic halos (Steidel et al. 1995; Mo 
& Miralda-Escude 1996) or in minihaloes formed prior to 
rcionization (Abel & Mo 1998). At the same time, the 
intensity of the cosmic UV background could be inferred 
from the observed brightness of the fluorescent emission. 

With present-day technology, the detection of fluo- 
rescent emission from high-redshift gas condensations is 
challenging, but not impossible. At z ~ 3, the intensity 
of the diffuse ionizing background (e.g. Haardt & Madau 
1996) corresponds to a Lya surface brightness of the or- 
der of 10 _20 erg cm -2 s _1 arcsec -2 . It is then not surpris- 
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ing that blind searches have only produced a number of 
null results (Lowenthal et al. 1990; Martinez-Gonzalez 
et al. 1995; Bunker, Marleau & Graham 1998). Positive 
fluctuations in the ionizing background can be used to 
increase the signal. For instance, clouds lying close to a 
bright quasar are exposed to a stronger UV flux (with 
respect to an "average" cloud) and are then expected to 
be brighter in fluorescent Lya. Very recently, Francis & 
Bland-Hawthorn (2004) presented a deep narrow-band 
search for Lya emission in a field which lies next to the 
quasar PKS 0424-131. Based on quasar-absorption-line 
statistics and on simple models for fluorescent emission 
(Gould & Weinberg 1996), they expected to detect more 
than 6 clouds but none were seen. These null results 
highlight the need for a more sophisticated analysis of 
fluorescent Lya emission in realistic environments. 

In this paper, we present accurate models of the fluo- 
rescent Lya emission from LLSs at redshift z ~ 3. Our 
study proceeds in three steps. First, we perform a hydro- 
dynamical simulation of structure formation to compute 
the cosmological distribution of the baryons at z = 3. 
A simple radiative transfer scheme is then used to prop- 
agate the ionizing radiation through the computational 
box and to compute the distribution of neutral hydro- 
gen and of recombinations. Finally, a three-dimensional 
Monte Carlo code is used to follow the transfer of Lya 
photons. As ionizing radiation, we first consider the dif- 
fuse background generated by the UV emission of galax- 
ies and quasars (Haardt & Madau 1996). We then dis- 
cuss an inhomogeneous case where the ionizing flux from 
a quasar (which lies in the foreground of the gas clouds) 
is superimposed to the uniform background. Our de- 
tailed numerical analysis shows that simplified models 
(e.g. Gould & Weinberg 1996) tend to overpredict the 
Lya flux emitted from optically thick regions. 

The structure of the paper is as follows. We describe 
our numerical techniques in §2 and present our results in 
§3 where we also discuss the implications of our analysis 
for present and future observations. Finally, we discuss 
the limitations of our approach in §4 and we conclude in 
§5- 

2. METHOD 

2.1. Cosmological simulation 

The formation and evolution of the large-scale struc- 
ture in a "concordance" ACDM cosmological model is 
followed by means of an Eulerian, grid based Total- 
Variation-Diminishing hydro+N-body code (Ryu et al. 
1993). We assume that the mass density parameter 
Oo = 0.3 (with a baryonic contribution Sib = 0.04), the 
vacuum-energy density parameter J1a = 1 — = 0.7 
and the present-day value of the Hubble constant con- 
stant H = 100 h km s" 1 Mpc" 1 with h = 0.67. The 
simulation is started at redshift z = 60 and follows 
the evolution of Gaussian density fluctuations character- 
ized by a primordial spectral index n = 1 and "cluster- 
normalization" as = 0.9 (with us the rms linear density 
fluctuation within a sphere with a comoving radius of 
8/i _1 Mpc). This is consistent with the most recent 
joint analyses of temperature anisotropies in the cos- 
mic microwave background and galaxy clustering (e.g. 
Tegmark et al. 2004 and references therein). We use a 
comoving computational box size of 10 /i -1 Mpc where 



the dark matter distribution is traced by 256 3 particles 
and the gas component is evolved on a comoving grid 
with 512 3 zones. The nominal spatial resolution for the 
gas (the mesh size) is ~ 20/i _1 kpc (comoving) with the 
mean baryonic mass in a cell being ~ 10 5 h^ 1 M Q . On 
the other hand, each dark matter particle has a mass of 

5 x 10 6 h~ x Mq. All the results presented in this work are 
derived from the z — 3 output of a simulation which does 
not include radiative cooling of the gas. The limitations 
of this assumption are briefly discussed in We defer 
a detailed analysis of the radiative case to future work. 

2.2. Radiative transfer of UV radiation 

In order to compute the distribution of neutral hydro- 
gen within a snapshot of the computational box, we need 
to simultaneously solve the radiative transfer problem for 
UV radiation and the rate equations describing the bal- 
ance between the ionization and recombination rates. 

For simplicity, we assume that hydrogen is in ioniza- 
tion equilibrium and use the "on the spot" approximation 
(Baker 1962): 

(l-a;)nH / dv-p^- dSlJ^Q) — xn H n c a B (T) (1) 

where hp, n c , x, nn, ov, T and «b respectively denote 
the Planck constant, the electron number density and 
the hydrogen ionized fraction, volume number density, 
ionization cross section, temperature and case B recom- 
bination coefficient (for which we use the fit by Hui & 
Gnedin 1997). The intensity of ionizing radiation per 
unit frequency and solid angle is given by 3 V (in erg 
cm~ 2 s _1 sr _1 Hz _1 ). The frequency integral in equa- 
tion Q extends from the hydrogen ionization threshold, 
/ipfo=13.6 eV, to a maximum frequency v up (which is, 
formally, infinite). A good approximation for our pur- 
poses is to assume v up = 4^o, (i.e. set the intensity 
of radiation to zero at frequencies above the ionization 
threshold for Hell). The motivation is twofold. First, 
nearly all the photons with v > 4 vq (which anyway con- 
tribute only a few per cent of the energy available for H 
ionization in the UV background) will be absorbed by 
He atoms (Haardt & Madau 1996). Second, Hell recom- 
bines faster than HI and the intensity of radiation at the 
Hell Lyman limit is typically lower than at i>q. Therefore, 
Hell is more easily shielded from the ionizing background 
with respect to HI (Miralda-Escude & Ostriker 1990). 
This implies that Hell-ionizing photons are absorbed in 
the outer regions of the gas concentrations where H is 
nearly fully ionized. In order to describe the hydrogen 
shielding layers we thus neglect Helll and assume that 
the neutral fraction of He coincides with 1 — x (Zheng 

6 Miralda-Escude 2002a). For a helium abundance of 
Y = 0.24, this corresponds to assuming n c = (3 x njj with 
(3 ~ 1.08 (see also fl2.3[) . Other than this, the presence of 
He atoms is neglected in equation {Q. Given that Hell 
recombinations produce Hi-ionizing photons and the rel- 
atively small number density of helium atoms and ions, 
this approximation should be reasonably accurate. Note 
that also recombination radiation from Helll can ionize 
HI. However, considered the different spatial distribution 
of Helll and HI discussed above and the characteristic 
Helll-recombination time scales, we neglect the small lo- 
cal corrections to the Hi-ionizing background deriving 
from this effect. 
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In each cell of the simulation, the diffuse ionizing back- 
ground is approximately described by following the ra- 
diative transfer along 6 "light-rays" which propagate par- 
allel (and antiparallel) to the main axes of the compu- 
tational box. With this numerical trick we can treat 
anisotropic backgrounds (created, for instance, by shad- 
owing effects) with a minimal request of CPU time (see 
Appendix lAl for a test of this approximation). Let us 
denote by Ti(v) the optical depth of a given cell along 
the i-th ray. This quantity is computed by integrating 
the product (1 — x) rva a v from a given starting location 
(a light source) in the box (see below) up to the first 
point of the cell crossed by the ray. The closest face of 
the cell is then exposed to a radiation field with intensity 
J^ n e" T *M , where J^" denotes the input ionizing radia- 
tion before it is filtered by the gas distribution in the 
box. Let us also indicate with At(v) = (1 — x) nu o v L 
(with L the cell size in physical units) the optical-depth 
variation within the cell measured along one of its prin- 
cipal axes. In order to implement a photon conserving 
scheme, we replace the left-hand side in equation Q with 
the quantity 

A~ 6 /-"up /in i _ -AtO) 

where the sum is taken over the six rays (labeled by the 
index i). This corresponds to the number of ionizing 
photons (per unit volume and time) which are deposited 
in a given cell by the six rays. To describe the diffuse 
UV background, we assume that jfy = J™ with J™ 
the intensity of radiation derived at z = 3 by Haardt & 
Madau (in preparation, hereafter HM) considering the 
emission from observed quasars and galaxies after it is 
filtered through the Lya forest 1 . We assume that under- 
dense cells are exposed to the full, isotropic background. 
On the other hand, overdense cells see an anisotropic ra- 
diation field which is computed by using equation (J2J) to 
propagate the input background starting from the sur- 
face where p = p. The intensity of radiation (and thus x) 
in each overdense cell depends on the ionized fraction of 
the surrounding region. To solve the non-local equations, 
we start our calculations by assuming that the whole sim- 
ulation box is optically thin (i.e. it is exposed to the in- 
put radiation field) and we iterate the radiative transfer 
and ionization-cquilibrium calculations until convergence 
(within 1%) is reached in each overdense cell. 

We use a similar approach to discuss the anisotropic 
radiation field generated by a quasar lying in the fore- 
ground of the simulation box along the observer's line 
of sight. For simplicity, we assume that the quasar lies 
distant enough from the simulated region that its emis- 
sion can be modeled as a train of plane waves imping- 
ing onto a face of the simulation box. We also assume 
that the quasar input spectrum is identical to that of 
the cosmic background. Given that J^ M is well de- 
scribed by a power-law of index -1.25 between Uq and 

1 This is obtained using the most recent results regarding the 
quasar luminosity function and cosmic evolution within a concor- 
dance cosmological model. It assumes that the galaxy escape frac- 
tion of ionizing raduation is / esc = 0.1 and that the energy spectral 
index for quasar radiation is a = 1.8. The resulting hydrogen ion- 
ization rate is a factor 1.16 smaller than in the models by Haardt 
& Madau (1996) used by Gould & Weinberg (1996). The spectrum 
is available at http://pitto.mib.infn.it/~haardt/refmodel.html. 



3^o, this is a sufficiently good approximation for our pur- 
poses (see also the extensive discussion in § I4.4f> . We 
then write the quasar ionizing flux (erg cm~ 2 s _1 Hz -1 ) 
as F v — nb J^ M 5u with <5y the Kronecker symbol and 
b a dimcnsionlcss constant. This is equivalent to using 
J"\ = 1.5 & J™ Su in equation I n this case, we 

compute the optical depth starting from the face of the 
simulation box which is first reached by quasar light (i.e. 
along the direction i = 1). 

A self-consistent calculation of the gas temperature re- 
quires a joint treatment of radiative transfer and hydro- 
dynamics which is still beyond present-day computing 
capabilities. Assuming that the photoionized gas is in 
thermal equilibrium, we find that T ~ 1 — 3 x 10 4 K for 
the typical densities in the shielding layers (100 < p/p < 
300). However, shock heating can easily drive the gas 
temperature to 10 5 ~ 7 K. This is particularly important 
for the low-density regions (p < 100 p) where cooling pro- 
cesses are inefficient and the shocked material remains 
hot (Theuns et al. 1998). In our analysis, we assume 
that T — 2 x 10 4 K everywhere. This is an excellent 
approximation for highly overdense regions (p > 100 p) 
where the cooling time is shorter than the Hubble time 
and the gas temperature rapidly approaches the equilib- 
rium solution (Theuns et al. 1998). Anyway, since the 
recombination coefficient «b has only a weak dependence 
on T, fixing the temperature to 2 x 10 4 K in the whole 
simulation box does not seriously affect our results. 

Note that, at T = 2 x 10 4 -ftT, the hydrogen recombi- 
nation timescale is t ICC — 2.26 (p/p) x 10 10 yr. Ioniza- 
tion equilibrium will approximately hold only where t Tec 
is shorter than the characteristic quasar lifetime (~ 10 8 
yr, Porciani, Magliocchetti & Norberg 2004), i.e. for 
p > 200 p. At lower densities, our assumption of ioniza- 
tion equilibrium will then overestimate the hydrogen ion- 
ized fraction. This is not a problem for our study since, in 
the vicinity of a quasar, the ionizing flux is strong enough 
to nearly completely ionize the low-density intergalactic 
medium. It is worth noticing, however, that regions with 
p < 200 p will emit their recombination radiation after 
the quasar has switched off and will not be detectable in 
a survey centred onto a bright quasar. 

2.3. The clumping factor 

Hydro-simulations have a finite spatial resolution and 
cannot describe the gas distribution on arbitrarily small 
scales. In other words, they provide a coarse grained rep- 
resentation of the density field. However, the hydrogen 
recombination rate scales proportionally to the square of 
the local (i.e. fine grained) number density and is sen- 
sitive to small-scale inhomogeneities (dumpiness) within 
a simulation cell. In order to keep track of this discrep- 
ancy, we re-write the mean recombination rate within a 
cell as 

x 2 pC(n H ) 2 a B (T) (3) 
where the average is taken over a simulation cell and 



denotes the clumping factor of the gas (we assume that 
different atomic species and ions have the same spatial 
distribution). In principal, the latter quantity can be 
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estimated by comparing simulations with different reso- 
lutions and consistent initial conditions. We assume that 
C is constant everywhere and we fix its value by imposing 
that the number density (per unit redshift) of LLSs in 
our simulation matches the observational data (Peroux 
et al. 2003). 2 This normalization procedure, which re- 
quires PC ~ 6, partially overcomes the limitations of our 
simulation (limited resolution and any missing physics). 
Note that the observational data constrain the product 
(3C so that there is no need to specify a priori the He 
ionization state as discussed in H2.2I 

2.4. Lya emission 

Using equation we compute the hydrogen recom- 
bination rate in each cell of the simulation. In order to 
convert this quantity into an emission rate for Lya pho- 
tons, we need to evaluate how many recombinations ul- 
timately lead to a 2 P — > 1 S transition. For T = 2 x 10 4 
K, nearly 44% of the atoms directly recombine to the 
ground level while 35% of the remaining cases ultimately 
produce excited atoms in the 2 S state which decays to 
1 S via two-photon emission (both fractions are weakly 
dependent on the gas temperature, see e.g. Osterbrock 
1989). Therefore, if the gas is optically thin to UV pho- 
tons, only a fraction e t hin = o^p/ckA ^0.36 (where a| p 
and a a denote the effective recombination coefficient to 
the 2P level and the case A total recombination coef- 
ficient, respectively) of the recombinations yield a Lya 
photon. However, in the optically thick case, continuum 
photons produced by recombinations to the ground level 
can be captured by neutral atoms and produce additional 
Lya radiation. The asymptotic yield in the extremely 
thick case (case B approximation, where no continuum 
photon can leave the cloud) is e t hick = «2p/«B ~ 0.65. 
We use this value to compute the emission rate of fluo- 
rescent Lya photons in the simulation box. 

2.5. Resolving the optical depth 

When we apply the method described above to our 
simulation, we find that the shielding layers (where the 
transition between optically thin and optically thick re- 
gions occurs) are poorly resolved (see Fig. Typically, 
they consist of very few cells which each correspond to 
an HI optical depth variation (at the Lyman limit) of 
At co h = At(z^o) > 1. However, for a proper treatment 
of the radiative transfer problem, more stringent require- 
ments on the grid spacing must be met. In particular, the 
Lya-emitting regions must be resolved with Ar co u < 1. If 
not, both the spatial distribution of recombinations and 
the escape probabilities of Lya photons along different 
directions (see H2.(jf) are spuriously altered. 

To solve this problem, we adaptively refine the Lya 
emitting regions by interpolating the original density and 
velocity fields of the input simulation. We use the so- 
lution of the radiative-transfer problem for the original 
(unrefined) grid to select the regions to interpolate and 
the factor of refinement. Given the memory limitations 
of the available machines, we use a 100 3 cells sub-box 
(which is particularly rich of structures) of the original 

2 Note that the spectral resolution of the observational data 
roughly corresponds to our box size. Therefore we can safely com- 
pute the hydrogen column density by integrating njji along the 
entire box. 




Fig. 1. — The fraction of all hydrogen recombinations happening 
in the simulation box which take place in cells whith a single-cell 
HI optical depth (at the Lyman Limit) Ar ce n > At; 011 is plotted 
as a function of AT; on . The function F rcc equals unity when Ar; on 
is equal to (or smaller than) the minimum optical depth of a single 
cell in the simulation box. Dashed and solid lines respectively refer 
to the original and the adaptively refined simulation boxes. Note 
that a proper treatment of the radiative transfer problem requires 
that the Lya photons are generated within cell with At cc u < 1 
(sec text). 



simulation and we interpolate every cell with a signifi- 
cant recombination rate (> 0.1% of the maximum) and 
Ar cc n > 1. The level of refinement is scaled proportion- 
ally to ATceii (up to a factor of 32 in each dimension) in 
order to have a subgrid of cells with Ar cc n < 1. Eventu- 
ally, we re-compute the radiative transfer for the adap- 
tively refined grid. Figure ^ shows that the fraction of 
recombinations originated in cells with Ar ce u > 1 de- 
creases from 30% to 7% as a result of this refinement. 
Moreover, in the finer grid, only a negligibly small num- 
ber of recombinations takes place in extremely thick cells 
(Ar ce n > 10) compared with 12% of the original grid. 

As discussed in H2.3I we account for unresolved sub- 
structure in our simulation box by using a non- vanishing 
clumping factor in the equation of ionization equilibrium. 
Density variations within a parent cell of the original sim- 
ulation due to the refinement procedure described above 
could, in principal, significantly contribute to the clump- 
ing factor. If this is the case, we should then adopt a 
value C < 6 for the refined simulation to reproduce the 
observed abundance of LLSs. We find that the clumping 
associated with the refinement is severe in the densest 
zones of the simulation (which typically lie in the self- 
shielded regions and do not contribute to the Lya flux) 
but amounts to only a few per cent in the most rapidly 
recombining cells. For these, we can then safely adopt 
C = 6 also for the refined box. 

Increasing the spatial resolution of the simulation com- 
plicates the radiative transfer of ionizing radiation gener- 
ated by recombinations. Equation Q assumes that every 
ionizing photon generated by a HII recombination is ab- 
sorbed in the same cell in which is generated. However, 
this is no longer a good approximation for the adaptively 
refined cells which are optically thin to UV radiation. In 



Fluorescent Lya emission 



5 



this case, ionizing photons generated by recombinations 
can be absorbed in a different cell with respect to where 
they are created. This process is too complicated to fol- 
low without an accurate radiative transfer scheme and 
we use equation JJJ also for the refined cells. How does 
this affect our results for the distribution of HI? First, the 
propagation of recombination radiation can slightly ex- 
tend (of a few cells) the thickness of the shielding layer 
of a gas cloud with respect to our results. The effect 
is probably more pronuciated in the outer shells where 
the gas density is lower. This should only redistribute 
the birth point of a small fraction of line photons. On 
the other hand, in the central part of the shielding layer 
(which contributes most recombinations) we expect that 
the flows of incoming and outcoming recombination ra- 
diation should nearly balance given that the hydrogen 
density shows little variations. In summary, our approxi- 
mated treatment of recombination radiation should only 
slightly modify the spectral energy distribution of the 
emerging Ly-a line 

2.6. Lya radiative transfer 

We now combine the results of the previous sections 
(namely, a set of arrays containing the Lya emission 
rate, the HI density and the gas velocity field as a func- 
tion of spatial position) to compute the spectra and the 
projected image on the plane of sky of the fluorescent 
sources. The radiative transfer of resonant Lya pho- 
tons is modeled using a three-dimensional Monte Carlo 
scheme analogous to that employed by Zheng & Miralda- 
Escude (2002b, see also Ahn, Lee & Lee 2001). The 
method follows a large number of photon trajectories as 
they are scattered within the HI density and velocity dis- 
tribution of the hydro-simulation. 

2.6.1. Emission of Lya photons 

We assume that Lya photons are isotropically emitted 
with frequency i>o in the frame of the recombining atoms 
(the natural linewidth is negligibly small for our pur- 
poses). In the cosmic frame (e.g. for an observer lying at 
the center of the simulation box and which participates 
to the free expansion of the universe), the frequencies 
of the resonant photons appear Doppler shifted by the 
projected velocities of the atoms along the photon tra- 
jectories. The velocity of a hydrogen atom with respect 
to the cosmic frame is given by the superposition of the 
Hubble flow with the bulk motion of the gas (i.e. the 
peculiar velocity of the fluid in the corresponding cell of 
the simulation) and a random thermal velocity: 

v = H(z)r + v gas + v th (5) 

with r the atom position with respect to the center of 
the simulation box. The component of v t h along the 
direction of the emitted photon is generated by extract- 
ing a Gaussian deviate out of a distribution with zero 
mean and dispersion er t h = (&b T/m^) 1 / 2 = 12.8 (T/2 x 
10 4 K) 1 / 2 kms -1 (with fee the Boltzmann constant and 
?71h the atomic mass). 

2.6.2. Absorption 

The photon frequency can be conveniently expressed 
in terms of the variable 



which measures the frequency shift from the Lya line 
center in units of the Doppler width, A = V^foCth/c, 
where c denotes the speed of light. The mean scattering 
cross section of Lya photons in the fluid frame is 

o-h ya (x) = Vtt fhya^- H(a,x) (7) 

where /Lya=0.416 is the Lya oscillator strength, r c = 
2.82 x 10~ 15 m is the classical electron radius and 

H(a,x) = - ^— -dy (8) 

Ti" J-oo (x-y) z +a z 

is the Hjerting-Voigt function. For the relatively low- 
densities we are interested in, atomic collisions are not 
important and the damping coefficient a can be ex- 
pressed in terms of the spontaneous decay rate T as 
a = r/(47rA) = 3.3 x 10~ 4 (T/2 x 10 4 K)~ 1 / 2 . 

We use equation (J7J| to determine the distance covered 
by each photon before it is scattered by an atom. We 
first extract a random deviate, R, from an exponential 
distribution function and then we integrate the product 
"■hi o~Lya {%) along the photon direction of motion until 
the resulting optical depth equals R. If the photon still 
lies within the computational volume, we select the ve- 
locity of the scatterer. Note that, in order to be able to 
absorb line radiation, an atom must have a velocity com- 
ponent along the trajectory of the incoming photon, vu , 
which closely matches the Doppler shift. From equation 
||SJ|, it follows that, in the fluid frame, xii — v\\/(^/2 er t h) 
is characterized by the following probability distribution 




We use the method presented by Zheng & Miralda- 
Escude (2002b) to generate deviates which follow this 
statistic. The perpendicular component of the ther- 
mal velocity in the scattering plane, x±, is then ex- 
tracted from a Gaussian distribution with a temperature- 
dependent dispersion as described above. 

2.6.3. Re-emission 

A new direction for the photon is then randomly se- 
lected according to a phase function, P(cos 9) (with 9 the 
scattering angle), determined by atomic physics. Res- 
onant scattering has an isotropic angular distribution, 
P = 1, while wing scattering is characterized by the 
Rayleigh phase function, P = 3(1 + cos 2 9) /A (Stcnflo 
1980). We find that the two angular distributions give 
consistent outputs. All the results presented in this work 
are obtained assuming isotropic re-emission. 

To determine the new photon frequency, we assume 
that the scattering process is coherent in the reference 
frame of the scatterer (partially coherent scattering). 
This is appropriate when the excited atom undergoes no 
collisions before re-emission and the radiative damping 
coefficient is small (Avery & House 1968). Both condi- 
tions apply to Lya radiation emitted by gas in the typi- 
cal conditions of the shielding regions in the intergalactic 
medium. Once the scattering angle and the photon ve- 
locity of the scatterer are specified, it is straightforward 
to compute the frequency shift of the re-emitted photon 
in the fluid frame: 

X = (Xin — %\\) + SET 1 1 COS-0 + X± SULIp (10) 
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Fig. 2. — Ly« spectrum emitted by a uniform slab with a mid- 
plane source with optical depth to. The results of our Monte Carlo 
code (solid histograms) are compared with an analytical approxi- 
mation (Neufeld 1990) which becomes exact in the limit to — > oo 
(dotted lines). A temperature of T = 10 K is assumed. 



where x m is the frequency shift of the incoming pho- 
ton and ip is the angle between the direction of the in- 
cident photon and the direction of the scattering atom. 
A Lorentz transformation is finally used to compute the 
frequency shift in the cosmic frame. 

The set of calculations described above is iterated until 
the photon escapes the computational box. 

2.6.4. Lya spectra 

To produce spectra (and broad-band images) of the flu- 
orescent emitters, we compute the surface-brightness of 
the computational box along the observer's line of sight 
(hereafter, the x-axis). At each scattering, the probabil- 
ity that a photon will be re-emitted along this direction 
is 

■^-P(cos6 x )e- T * (11) 

where 9 X is the angle between the incoming photon and 
the x-axis and t x denotes the Lya optical depth of the 
scattering site along the observer's line of sight. 3 For 
each photon and for each scattering, we sum this quan- 
tity to a counter in correspondence of the projected po- 
sition of the scattering site and of the photon frequency. 
We thus obtain a three-dimensional array containing the 
surface brightness of fluorescent Lya photons as a func- 
tion of 2 spatial coordinates plus frequency. Note that 
a simulated photon tends to remain for many scatter- 
ings in a rather small region before it eventually escapes. 
This means that photons contribute only to a few pixels 
surrounding their emission site. 

Following Zheng & Miralda-Escude (2002b), we test 
our implementation of the Monte Carlo scheme against 
the analytical approximation by Neufeld (1990) for the 
optically thick, plane-parallel case. Figure [21 shows that 
our code accurately reproduces the analytical solution 

3 This optical depth includes the effects of neutral hydrogen lying 
in the foreground of the computational box. 



which becomes exact in the limit of extremely large op- 
tical depths. 

3. RESULTS 

In order to have an acceptable compromise between 
spectral resolution and CPU time, we only apply the 
Monte Carlo radiative transfer to the adaptively refined 
grid corresponding to a 100 3 region of the original simula- 
tion box. To achieve a good signal-to-noise ratio, we gen- 
erate 10 6 photon trajectories for every simulation. We 
thus obtain high resolution spectra for each pixel of the 
resulting image that can be combined to simulate slit, 
line-emission integral field or broad-band observations. 

3.1. Diffuse background and static gas 

We first discuss the ideal case of a static gas distribu- 
tion illuminated with a uniform and isotropic backround 
of ionizing radiation. This is obtained by artificially set- 
ting to zero the velocity field of the gas within our refined 
box. 

In the left panels of Figures and we respectively 
show the HI column density distribution and the broad- 
band images (~ 90 A in the oberved frame, centered at 
A = (1 + z) ■ 1216 A = 4864 A) of the selected region 
illuminated with the diffuse UV background. The color 
code in Figure 01 gives the fluorescent Lya emission rate 
(photons per unit time, surface and solid angle) in units 
of the impinging rate of ionizing photons times e (i.e. the 
fraction of the recombinations yielding a Lya photon): 

riva jHM 

Rum = ethick / -r—fo = 2.44 x 10 4 cm" 2 s" 1 sr" 1 

with e t hick — 0.65. For an observer at redshift z = 0, 
this corresponds to a Lya surface brightness of 

SBhm = 3.67 x 10~ 20 erg cm" 2 s" 1 arcsec" 2 . (13) 

The brightest fluorescent sources correspond to compact 
gas clouds with a meatball topology. This is because 
the diffuse UV background is bright enough to fully ion- 
ize gas concentrations with p < 100 p. In general, the 
shielding regions either lie within virialized structures or 
correspond to dense gas shells which are accreting onto 
collapsed objects. As we will see below, the velocity field 
of the infalling gas produces specific signatures in the 
Lya spectra. 

The compact fluorescent sources lie along the filaments 
and sheets which characterize the distribution of neutral 
hydrogen on cosmological scales. For ease of reading, 
we label the three largest structures (which each have 
a diameter of ~ 0.4 comoving Mpc) with the letters A, 
B and C (see Fig. 0J. Cloud C is composed of two 
sub-units and is a part of an elongated structure which 
extends towards cloud B. Similarly, a filamentary plume 
of gas bridges clouds A and B. 

4 There is some observational evidence that the UV background 
at z = 3 is dominated by quasar emission with a negligible con- 
tribution from star-forming galaxies (e.g. Scott et al. 2000). In 
this case, the models by Haardt & Madau (in preparation) give 
Rum = 1-88 x 10 4 cm -2 s _1 sr _1 . The spectral shape of the UV 
background between and 4i/o is nearly identical to the general 
case discussed in the main text. Therefore, our predictions for 
the surface brightness of fluorescent sources can be simply scaled 
down by 30% if future observations will prove that galaxies do not 
significantly contribute to the ionizing background at z = 3. 
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Fig. 3. — Column-density distribution of neutral hydrogen at z = 3. In the left panel, the gas is exposed to a diffuse UV background 
generated by the population of galaxies and quasars. In the right panel, the ion izin g flux from a foreground quasar, located a short distance 
in front of the region and corresponding to a boost factor 6 = 6 (see equation 1141 ) is superimposed to the diffuse background. 
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Fig. 4. — Broad band images (~ 90 A, centered on 4864 A) of fluorescent Lya emission at z = 3 for static gas clouds (left) and accounting 
for the gas velocity field (right). Both images correspond to the colum density distribution presented in the left panel of Figure The 
white boxes indicate the location of the slit spectrographs used to obtain the energy distributions presented in Figure [7| 



Simple reasoning based on the plane-parallel model 
for line transfer suggests that, in the absence of photon 
sinks (e.g. dust), self-shielded (isotropically-illuminated) 
objects should shine with a surface brightness of SBhm 
(Hogan & Weymann 1987; Gould & Weinberg 1996). In 
our static simulation (Fig. 0] left panel), the SB of self- 
shielded objects closely matches the predictions of this 
simple plane-parallel model. The SB distribution in the 
simulation (dotted histogram in Fig. EJ) shows a nar- 
row peak at this expected value. In general, the SB 

1/2 

scales proportionally to iV HI in the optically thin re- 
gions and asymptotically approaches its maximum value 
for self-shielded objects (see the top-left panel in Fig. 
Hj}. The brightest lines of sight in fluorescent Lya corre- 



spond to optically thick systems with column densities 
Nni > 10 18 cm 2 which are thus associated with LLSs and 
DLAs. All the photons of the ionizing background are 
converted into Lya radiation within the shielding layers 
of these optically thick systems. In the absence of other 
sources of ionizing radiation, it is impossible to produce 
a stronger Lya flux. This explains why the brightest ob- 
jects in the left panel of Figure 01 have a uniform SB and 
sharp boundaries which correspond to the regions with 
N m ~ 10 18 cm 2 in the left panel of Figure 

3.2. Diffuse background and realistic gas velocities 

We are now ready to discuss the more realistic case 
where we include the gas velocity field of the hydro- 
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Fig. 5. — Surface brightness distribution of fluorescent sources 
ionized by a diffuse UV background (solid), by the additional con- 
tribution of a quasar with "boost" factor 6 = 2 (dashed) and 6 = 6 
(dot dashed). The dotted line is analogous to the solid one but is 
obtained by artificially setting to zero the gas velocity field. 



simulation. The corresponding Lya emission rate is 
shown in the right panel of Figure 01 The overall pattern 
is similar to the static case, but a number of striking 
differences are noticeable. Namely: i) the SB of self- 
shielded objects is no longer uniform (e.g. the right-hand 
side of Cloud A is nearly a factor of 2 fainter than the 
left-hand side); ii) the boundaries of the emitting regions 
are less sharp and self-shielded objects are surrounded by 
large, low-SB halos; Hi) self-shielded objects can be sig- 
nificantly fainter (or, very rarely, brighter) than in the 
static case. 

The top- right panel in Figure El shows that the gas ve- 
locity field introduces additional scatter into the SB - 
A'hi relation with respect to the static case. The bright- 
est lines of sight still correspond to Nm > 10 18 cm 2 but 
now two regions with the same column density can be as- 
sociated with brightnesses which differ up to a factor of 5. 
In consequence, the SB distribution of optically thick re- 
gions is broader and it is slightly shifted to fainter fluxes 
with respect to the static case (see the peak of the solid 
histogram in Fig. |SJ. We find that the median SB of 
the self-shielded objects amounts to nearly 75% of the 
value predicted by Gould & Weinberg (1996). At the 
same time, a larger fraction of the sky has SB < SBhm 
compared to the static case, (the power-law part of the 
solid histogram in Fig. As we will show below, this 
excess is caused by foreground scattering of the Lya pho- 
tons and is related to the presence of extended Lya halos 
around self-shielded objects. 

A better understanding of the "velocity-field effect" 
can be achieved by comparing the spectra of the fluo- 
rescent emission in the static and in the general case. 
In the left and central panels of Figure [7| we show the 
corresponding spectral energy distributions of the Lya 
photons. These have been obtained positioning four slit 
spectrographs (width ~ 0.9 arcsec and variable length) 
on top of the three brightest sources as shown in the right 
panel of Figure 0] In a static gas distribution, spectra 




l°g(N HI /cm- 2 ) log(N„,/cm- 2 ) 

Fig. 6. — The Lyo surface brightness of each pixel of the simu- 
lated images is plotted against the corresponding column density 
of neutral hydrogen. In the top panels, the intergalactic medium is 
ionized by a diffuse background. In particular, the top-left frame 
refers to a static gas distribution. In the bottom panels, a quasar 
with boost factor 6 = 2 (bottom-left) and 6 = 6 (bottom-right) 
is superimposed to the diffuse background. Dotted lines mark the 
expected SB for a plane-parallel slab while dashed lines indicate 
the minimum column density for LLSs (short-dashed) and DLAs 
(long-dashed). 



have a characteristic double humped shape and are sym- 
metric with respect to the line center. On the other hand, 
in the general case the energy distribution is no longer 
symmetric. In fact, particular configurations of the ve- 
locity and density fields are able to strongly suppress one 
of the wings of the Lya line and significantly lower the 
observed SB of the self-shielded objects. In the particular 
case of Cloud A, a low-density concentration of neutral 
hydrogen is infalling onto the Lya emitting region. The 
relative velocity (along the line of sight) corresponds to 
~ 4<7th and thus to a very high optical depth. Therefore, 
most of the photons that, in the static case, leave the 
shielding layers along the line of sight in the red Doppler 
wing will scatter within the infalling cloud and escape in 
other directions loweing the observed surface brightness. 
These photons will then form the extended Lya halos 
which surround the brightest objects in Figure The 
phase-space distribution of neutral gas in the vicinity of 
the emitting regions thus plays a fundamental role in re- 
shaping the Lya spectral energy distribution. In broad 
terms, infalling material diminishes the red wing of the 
spectrum, while gas which is receding from the emitting 
region (which could also mean that the shielding layer is 
infalling onto a central object more rapidly than the sur- 
rounding gas) damps out the blue peak of the spectrum. 
We find that the velocity dispersion of the gas within 
the regions crossed by the Lya photons broadens the red 
and blue peaks of the spectrum by up to 10 A. On the 
other hand, when both peaks are detectable, their sepa- 
ration is nearly independent from the detailed properties 
of the emitting regions and is set by the thermal veloc- 
ity dispersion of the original cloud. In fact, in analogy 
with the plane-parallel case, the escape probability of a 
Lya photon peaks at a frequency which only depends on 
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Fig. 7. — Two-dimensional spectra obtained "observing" our simulations with the slits shown in the right panel of Figure^H Different 
columns refer to different simulations. In p articular, from left to right: diffuse backgr ound and static gas distribution f ij3.ll . diffuse 
background and realistic gas velocities (Jl^l, quasar with 6 = 6 plus diffuse background f H3.3I , The labels in the top left corner indicate 
a particular slit spectrograph as in FigureBI 



the optical depth of its emission site and on the temper- 
ature of the medium. For a typical self-shielded cloud 
(tll > 1), the spectrum peaks at ~ ±4x which, for the 
assumed temperature, corresponds to a separation of ~ 8 

A. 

The two-dimensional spectra shown in the central 
panel of Figure [7| clearly show that the gas velocity field 
within and in the vicinity of the shielding layers has a 
complicated structure which does not show the charac- 
teristic pattern of ordered rotation or symmetric infall 
considered by Zheng & Miralda-Escude (2002b). 

3.3. Quasar plus diffuse background 

We now discuss a case of anisotropic illumination, ob- 
tained by superimposing the ionizing flux from a quasar 
to the diffuse UV background. The quasar is imagined 
to lie a short distance in front of the computational box 
as seen by us, and thus enhances the UV illumination 
experienced by faces of gas clouds exposed to it. Note 
that the "boost" factor b (defined in H2.2JI is determined 
by the intrinsic luminosity of the quasar and by its ac- 
tual separation from the simulated region. The definition 
above can be generalized to any given quasar spectrum 
using the emitted rate of ionizing photons. At a physical 
distance r from a quasar with monochromatic luminosity 



b = 15.2- 



0.7 



10 30 erg s" 1 Hz" 



IMpc 



(14) 



The resulting Nm distribution (assuming a boost factor 
b = 6) is shown in the right panel of Figure |21 The 
corresponding broad-band image (obtained accounting 
for gas velocities) is presented in the left panel of Fig- 
ure |H1 As expected, the self-shielded regions (and thus 
the fluorescent sources) are smaller with respect to the 
isotropic background case due to the extra-ionizing radi- 
ation produced by the quasar. This also makes the fluo- 
rescent sources brighter (dot-dashed histogram in Fig. [SJ 
since more recombinations will be produced to balance 
a stronger ionization rate. Based on the (plane-parallel) 
slab model, where Lya photons are emitted following a 
cosine law (Gould & Weinberg 1996), one would have 
naively expected an increase in the Lya surface bright- 
ness towards the observer by a factor 1 + 6 = 7 with re- 
spect to the diffuse background case. 5 However, Figures 
[S] and clearly indicate that the slab model overstimates 

5 This holds for normal incidence. In general, the surface bright- 
ness of a slab which forms an angle 9 with the incident quasar flux 
corresponds to a factor 1 + b cos 9. 
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Fig. 8. — Broad band images (~ 90 A, centered on 4864 A) of fluorescent Lya emission from our simulation box at z = 3. In both 
cases, the intergalactic medium is ionized by a diffuse background and by a quasar with boost factor 6 = 6. The image on the left (right) 
is obtained by observing our simulation box from a line of sight parallel (perpendicular) to the direction of quasar illumination. The left 
frame corresponds to the column density distribution presented in the right panel of Figure l3l 



the SB of the self-shielded objects. This is not due to 
shadowing effects. In fact, the attenuation of the quasar 
flux by diffuse gas lying in front of the fluorescent clouds 
is generally negligible. Comparing with a static simu- 
lation, we also find that gas motions can only explain 
a small part of this discrepancy. In fact, in the pres- 
ence of a quasar, foreground scattering is reduced due 
to the lower neutral fraction present in low density gas 
and broad-band images tend to be more uniform than 
in the case of isotropic illumination. On the other hand, 
the slab approximation no longer applies when the size 
of the shielding layers is comparable with the radius of a 
cloud. In this case, Lya photons produced at a particular 
point leave the cloud with a different angular distribu- 
tion with respect to the plane-parallel case. For approx- 
imately spherical clouds and in the presence of uniform 
illumination, this effect is suppressed for symmetry rea- 
sons. However, when the ionizing flux is anisotropic, the 
Lya SB does depend quite strongly on the geometry of 
the shielding layers. 

To study how the SB of self-shielded objects along the 
quasar direction, SBm, depends on the impinging flux, we 
performed a series of simulations with increasing b. Our 
results are summarized in Figure [!3 where we express 
SBy in terms of an "effective boost factor" defined by 

SB,, = (1 + b eG ) SBhm ■ (15) 

Points with errorbars mark the 25 th , 50 th and 75 th per- 
centiles of l+6 c ff among the DLAs. The solid line repre- 
sents the best-fitting relation 

1 + b cff = 0.74 + 0.50 6 89 , (16) 

while the dashed line shows the predictions of the slab 
model. Note that the geometric effect becomes more and 
more important as b is increased. 

Where do the "missing" Lya photons go? In the right 
panel of Figure|Sl we show the fluorescent emission along 



a line of sight perpendicular to the direction of quasar il- 
lumination (assuming b = 6 as in the left panel). In 
this case, the plane-parallel model predicts that the self- 
shielded objects should emit at SBhm- In our simula- 
tions, however, the shielding layer deeply penetrates in 
the clouds along the quasar direction and the slab model 
does not apply. In consequence, self-shielded objects 
are much brighter than a slab along this line of sight. 
Typically, SB ± ~ 0.5 SB), for b > 1 while SB ± ~ SB,, 
for b <C 1. In other words, Lya photons generated by 
the quasar ionizing flux are emitted within a wide solid 
angle. As a consequence of this partial isotropization, 
self-shielded clouds are fainter than expected (based on 
the slab approximation) along the quasar direction and 
brighter in the perpendicular directions. 

Finally, in the bottom panels of Figure El we show the 
SB - ATjji scatterplot for anisotropic illumination (when 
observer, quasar and the simulation box are aligned). It 
is worth noticing that, while the SB keeps nearly con- 
stant for LLSs, on average, it steadily increases with TVhi 
for DLAs. This phenomenon can be explained by as fol- 
lows. Let us assume that self-shielded objects are nearly 
spherically symmetric. Then, i) the ionizing flux from 
the quasar depends on the incident angle with respect 
to the local density gradient in the clouds; ii) this co- 
sine approaches unity for the central projected regions of 
self-shielded objects; iii) the column density reaches the 
highest values along these lines of sight. 

3.4. Size distribution of Lya sources 

Knowing the size distribution of fluorescent Lya 
sources is fundamental to planning an observational cam- 
paign for their detection. Regrettably, our refined box 
is too small (its size being ~ 2/i -1 comoving Mpc) to 
provide a statistically representative sample of optically 
thick sources. On the other hand, performing the line 
transfer on the 10 ft. -1 Mpc box would require an exces- 
sive amount of computer time. For these reasons, we 
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Fig. 9. — Lya surface brightness of optically thick clouds [ex- 
pressed in terms of the effective boost factor defined in equation 
1151 1 as a function of the impinging quasar ionizing flux [expr esse d 
in terms of the quasar boost factor, b, defined in equation 1141 1. 
Points with errorbars denote the 25 th , 50 th and 75 th percentiles of 
b e ff for the DLAs. T he so lid line represents the best-fitting rela- 
tion given in equation 1161 . Predictions of the static, plane-parallel 
model are plotted with a dashed line. 



decided to propagate only the ionizing radiation through 
the 10 h~ x Mpc box and to use the scatterplots in Fig- 
ure [S] to convert the neutral-hydrogen column densities 
into Lya fluxes. In fact, independently of the value of 
b, all lines of sight with TVhi > 10 18 cm -2 are approx- 
imately associated with a constant Lya surface bright- 
ness (within a factor of 2 uncertainty caused by the gas 
motion and cosine effects discussed above). We then 
adopt this threshold value to derive the size distribution 
of fluorescent objects. In Figure we present our re- 
sults for an isotropic ionizing background (b = 0). Solid 
and dashed histograms respectively refer to objects with 
Nm > 10 18 cm -2 and to DLAs. It is worth remembering 
that we fixed the value of the clumping factor in our simu- 
lation so as to reproduce the observed sky covering factor 
of LLSs. In consequence, if a significant fraction of the 
real systems have a characteristic size which is smaller 
than our numerical resolution, our simulation will over- 
predict the number of large systems in order to preserve 
the required normalization. 

In Figure 1111 we plot the number density of self- 
shielded objects as a function of b. We use three different 
thresholds for the source size: 3 (which corresponds to 
barely resolved objects), 20 and 80 arcsec 2 . In all cases, 
the number of sources rapidly drops with increasing b. 
In fact, higher values of b characterize regions which are 
closer to a given quasar (see eq. H140 ) and, obviously, 
correspond to a lower number density of self-shielded ob- 
jects. 

From this figure, it is also possible to determine the 
number density of sources which are brighter than a cer- 
tain threshold value (1 + 6 m in)SBHM- Let us consider 
a Lya source which is optically thick to ionizing radi- 
ation at a given distance from a quasar. Let us also 
imagine that we can move the cloud towards the quasar 




A (arcsec 2 ) 

Fig. 10. — Differential size distribution of optically thick clouds 
(solid histogram) and DLAs (shaded histogram) in a simulation 
where the intergalactic medium is exposed to a diffuse ionizing 
background. The shaded histogram has been slightly shifted in the 
horizontal direction to improve readibility. 




1 + b 

Fig. 11. — Physical number density of fluorescent Lye? sources 
(with size A indicated by the labels) as a function of the impinging 
quasar fl ux [ expressed in terms of the boost factor b defined in 
equation 1141 1. Errorbars are derived assuming Poisson statistics. 
Triangles and circles have been slightly displaced in the horizontal 
direction to improve readibility. 



thus increasing the b factor. As long as the cloud keeps 
optically thick, SB|| monotonically increases. However, 
there exists a particular value of the boost factor, 6 SS , at 
which the cloud is no longer able to self-shield. There- 
fore, for b > b ss , SB || keeps roughly constant. 6 Thus, 
the number of self-shielded objects at a given b coincides 
with the number of sources (which are not necessarily 
optically thick) with SB > [1 + b e s(b)] SBhm- In other 

6 The fraction of recombinations yielding a Lya photon decreases 
from ethick ~ 0.65 to e t hin ~ 0.36 when a cloud becomes optically 
thin. Therefore, we expect a fully ionized cloud to be a factor of 
~ 2 fainter in Lya with respect to the optically thick case. 
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words, the number of sources which are brighter than a 
given threshold can be computed with the following pro- 
cedure. First, convert the threshold SB into an effective 
boost factor, 6*^ r . Second, invert equation (|Ttj)l to find 
the value of 6 m i n such that b e g(b m i n ) = b^. Third, use 
b = b mm in Figure HTl to determine the number density of 
the sources. Fourth, use equation l|14fl to find the volume 
within which it is possible to have b > 6 m m- 

The variation of the number density of fluorescent 
sources as a function of b is somehow related to the prox- 
imity effect. Hydrogen clouds in the vicinity of a quasar 
are strongly ionized and emit fluorescent radiation. In 
other words, the missing absorption systems which de- 
termine the proximity effect can be detected in emis- 
sion through their recombination radiation. Therefore, 
in analogy with studies of the proximity effect (e.g. Ba- 
jtlik, Duncan & Ostriker 1988; Scott et al. 2000), the 
number density variation of fluorescent sources around 
a quasar can be used to infer the intensity of the UV 
background. In this case, reliable models of fluorescent 
emission are fundamental to convert the observed counts 
into a background amplitude. 

3.5. Comparison with recent observational data 

We can use the above to compare the predictions 
of our models with the recent observational results 
by Francis & Bland-Hawthorn (2004, hereafter FBH). 
These authors performed a deep narrow-band search 
for fluorescent Lya emission in the vicinity of the 
z = 2.168 quasar PKS 0424-131 (X LL = 1.67 x 
10 30 erg s _1 Hz _1 , a ~ 0.7). At the 5c confidence 
level (corresponding to a surface brightness of 4.7 x 
10~ 19 erg cm -2 s _1 arcsec -2 for sources larger than 
lOOarcsec 2 and to 9.6 x 10 _18 erg cm -2 for unresolved 
sources) no source was been detected. Based on the ob- 
served abundance of LLSs, FBH expected to find ~ 6 
fluorescent clouds with a size of 100 arcsec 2 . This esti- 
mate, however, does not take into account the ionizing 
radiation from the quasar. 

Assuming that our results at z — 3 are approximately 
valid at the quasar redshift, 7 we find that the sensitiv- 
ity limits of FBH correspond to b m - m ~ 11.2 for sources 
which are larger than 100 arcsec 2 . Assuming that the 
ionizing backround keeps roughly constant in the red- 
shift interval 2 < z < 3, from equation l|14|l we find 
that this corresponds to a distance from the quasar of 
'"max ~ 1.5 (physical) Mpc. This is the maximum dis- 
tance from the quasar at which an optically thick cloud 
could have been detected. Based on our simulations, we 
expect to find, on average, 1 — 2 objects within a sphere 
of radius r max around the quasar. However, FBH limited 
their search to distances smaller than 1 Mpc thus reduc- 
ing the sampled volume by a factor 3.4 with respect to 
the theoretical limit. In this case, the expected num- 
ber of sources ranges between 0.3 and 0.6. Therefore, 
the probability of detecting (at least) one source with 
one single observational run is 0.25 < P < 0.45 (assum- 
ing Poisson statistics) . Our simulations clearly show that 
the results by FBH are thus perfectly consistent with our 
understanding of the intergalactic medium at high-z. 8 

7 We simply assume that the Lyct surface brightness scales as 
(1 + z)- 4 , i.e. SBO = 2.168) = 2.54 SB(z = 3). 

8 The detection limit for unresolved sources is less interesting. 



There are caveats to the simple calculations descrived 
above. For instance, we have assumed that all the Lya 
sources lying within a distance r m i n from the quasar are 
brighter than & m ; n . This assumption tends to overesti- 
mate the number of detectable objects. In fact: i) self- 
shielded clouds lying in front of the quasar are much 
fainter and hardly detectable (their SB actually depends 
on the angle between the line of sight and the quasar di- 
rection); ii) fully ionized clouds in the foreground of the 
quasar tend to be a factor of ethick/ethin — 1-8 fainter 
than assumed above; iii) if the age of the quasar is 
shorter than the hydrogen recombination timescale no 
fluorescent sources will be detectable (see the discussion 
at the end of §2.2). On the other hand, we have as- 
sumed that our simulation box is representative of the 
gas distribution surrounding a quasar. Since optically 
selected quasars at high-z tend to sit within the most 
massive dark-matter halos formed at that epoch (Por- 
ciani et al. 2004), it is reasonable to expect that mat- 
ter clusters (and moves with larger peculiar velocities) 
around them. Therefore, we could have underestimated 
the number densities of fluorescent sources lying close to 
a quasar. 

Despite of the approximations listed above, we believe 
that our results provide the most reliable estimates for 
the abundance of fluorescent Lya sources at high-z car- 
ried out to date. Optimized sampling strategies are cer- 
tainly required to observe these objects. Our simulations 
then represent a fundamental tool to plan observations 
around a given quasar. 

4. UNCERTAINTIES 

4.1. Radiative cooling 

While numerical simulations are a useful tool to guide 
our understanding, they cannot be considered a perfect 
model of reality. A potential limit of our simulation is 
the lack of radiative cooling. While this is not a con- 
cern for the diffuse intergalactic medium, it becomes wor- 
risome for highly overdense regions. Fluorescent Lya 
sources have intermediate overdensities [fijp ~ 200) and 
are likely to be in equilibrium with the UV ambient radi- 
ation. We thus expect them to be mildly affected by cool- 
ing processes. In any case, most of the results discussed 
in this paper are nearly independent of the details of the 
gas distribution. Both the velocity and the geometric ef- 
fects will be present anyway. On the other hand, the size 
distribution of the sources might be more affected by the 
cooling processes. It is also worth stressing that there is 
no way of accounting for the effects of cooling and heat- 
ing in a self consistent way. In fact, simulations where 
radiative transfer is fully coupled with hydrodynamics 
are still not viable with current supercomputers. More- 
over, other poorly understood processes (like energy and 
momentum feedback) play an important role. Therefore, 
it is hard to judge the level of approximation that current 
simulations including gas cooling may reach. 

4.2. Resolution and sub-structure 

Is the finite resolution of our simulation affecting our 
results? Multiple metal lines are often associated with 

It corresponds to fe m i n ~ 388 (assuming that equation 1161 can 
be extrapolated to such high values of b) and r max ~ 250 kpc. 
The number of expected sources in the associated volume of ~ 
0.06 Mpc 3 is therefore negligibly small. 
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single DLAs (e.g. Prochaska & Wolfe 1997) thus sug- 
gesting the presence of a clumpy medium. Unresolved 
sub-structure in our simulation might reduce the velocity 
effect and modify the outcoming spectra. The adopted 
value of C implies that at least 1/6 of the volume is in 
dense clumps. If these sub-structures have a diameter 
which is comparable to the cell size, we only expect a 
minor modification of our results. In fact, the gas ve- 
locity in the simulation should closely approximate the 
motion of these clumps. The only effect is then a slight 
Doppler-shift of the entire Lya spectrum of each cell. 
The opposite case, where each cell contains a large num- 
ber of small substructures (Abel & Mo 1998), can be ap- 
proximately discussed by considering an additional con- 
tribution to the thermal velocity dispersion. Assuming 
a value of a t h ~ 50kms _1 (corresponding to roughly 
half the virial velocity of the host halos, Haehnelt, Stein- 
metz & Rauch 1998), we find that the velocity effect 
is be strongly suppressed and the spectra are in better 
agreement with the slab model. Note, however, that the 
existence of a sea of small subclumps is disfavored by 
observations. In fact, this scenario would produce broad 
absorption features instead of the multiple metal systems 
associated with single DLAs (cf. Haehnelt, Steinmetz & 
Rauch 1998; McDonald & Miralda-Escude 1999). In any 
case, the velocity dispersion within our simulated DLAs 
is of the order of 100 km s -1 , in good agreement with 
observational data. 

4.3. Additional sources and dust 

Beyond fluorescent emission, additional Lya radiation 
might be produced in the inner regions of the clouds. 
Within gravitationally collapsed objects, the gas tends 
to dissipate its internal energy by emitting line photons 
(e.g. Haiman, Spaans & Quataert 2000; Fardal et al. 
2001; Furlanetto et al. 2005). Similarly, internal star 
formation could act as a copious source of Lya photons, 
but, whereas fluorescent emission is expected to extend 
over several tens of kiloparsec, the Lya emission from 
star forming region should be more concentrated near 
the centers of galaxies (Furlanetto et al. 2005). We have 
focused here on the fluorescent emission generated by 
recombinations and these extra sources of line photons 
are not considered in our analysis. We will present a 
comprehensive model of Lya emitters in a future paper. 

At the same time we did not consider the destruction 
of Lya photons by dust grains. Little is known about 
the dust properties within the intergalactic medium at 
z ~ 3 even though there is some evidence for the presence 
of dust in DLAs (Fall & Pei 1993). Nevertheless, the 
associated absorption of fluorescent photons is likely to 
be minimal due to the relatively low iVni of the shielding 
layer (e.g. Gould & Weinberg 1996). On the other hand, 
absorption is expected to be more severe for Lya photons 
produced close to and within the star-forming regions 
where dust is likely to be more abundant and the Lya 
escape probability is lower. 

4.4. The UV background 

Our results are based on the Haardt-Madau model for 
the UV background. At z = 3, the adopted model has an 
amplitude at the Lyman limit of J(vo) = 4.04 x 10 -22 erg 
cm~ 2 s _1 sr Hz and corresponds to an hydrogen ion- 
ization rate T = 1.15x 10~ 12 s _1 . Even though on the low 



side, this is consistent with recent observational studies of 
the proximity effect, which give J(vq) = 7.0^%% X 10~ 22 
erg cm~ 2 s" 1 sr^Hz- 1 and T = 1.9±{% x 10~ 12 s" 1 
(Scott et al. 2000). Our results should then considered 
conservative, because they account for a UV background 
intensity which is a factor of 1.75 lower than current ob- 
servational estimates. Note, however, that a different 
normalization of the UV background cannot change our 
predictions for the number density of fluorescent sources. 
In fact, we determined the clumping factor of the gas by 
imposing that the mean number of LLS in the simulation 
matches the observational data. Therefore, a stronger 
ionizing radiation field would require a larger clumping 
factor to fit the observed counts. If both J^ M and C 
are multiplied by the same multiplicative factor, Aj, the 
solution of the ionization equilibrium does not change. 
In fact, this is equivalent to multiply both sides of equa- 
tion UJ by a constant factor. Therefore our results for 
the number density and size distribution of fluorescent 
sources are robust with respect to the amplitude of the 
UV background. Note, however, that the Lya surface 
brightness of optically thick cloud would be higher by a 
factor A j with respect to equation H13fl . Consistently, 
the boost factor in equation (|14|l would decrease by a 
factor Aj while equation (|16(l would remain unaltered. 

Some (most likely minor) modification of our number 
counts is instead expected if the spectral energy distribu- 
tion of the UV background significantly differ from the 
assumed one. This crucially depends on the relative con- 
tribution of galaxies and quasars to the cosmic ionizing 
background. Similarly, the spectral energy distribution 
of quasar radiation plays an important role. For simplic- 
ity, we assumed that the quasar spectrum is identical to 
that of the cosmic background. This corresponds to a 
energy index a ~ 1.25 i.e. to a much flatter spectrum 
than inferred from quasar observations (1.5 < a < 1.8). 
What could be the effect of this assumption? The aver- 
age mean free path of ionizing photons from a spectrum 
with a = 1.8 is only a factor of 20% smaller than for 
a = 1.25. Therefore, radiation from a steep quasar spec- 
trum should penetrate a bit less within hydrogen clouds 
with respect than our models. This might slightly modify 
the emerging Lya spectra and reduce the importance of 
the geometric effect. Note, however, that quasar radia- 
tion will be partially filtered by the IGM before reaching 
the fluorescent clouds. Thus, a spectrum with a — 1.8 at 
emission will be transformed into a flatter distribution in 
the shielding layer of a cloud. 

5. SUMMARY 

We have presented a new method to produce realis- 
tic simulations of fluorescent Lya sources at high red- 
shift. We started by simulating the formation of baryonic 
large-scale structure in the ACDM cosmology. A sim- 
ple radiative transfer scheme was then used to propagate 
ionizing radiation through the computational box and to 
derive the distribution of neutral hydrogen. Finally, the 
transport of Lya photons generated by hydrogen recom- 
binations was followed using a three-dimensional Monte 
Carlo code. As ionizing radiation, we first considered the 
smooth background generated by galaxies and quasars. 
Then, as a second case, we superimposed to the back- 
ground the UV flux produced by a quasar lying in the 
vicinity of the simulation box. 
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Our detailed numerical treatment improves upon pre- 
vious work which was either based on rather crude ap- 
proximations for the transfer of resonantly scattered ra- 
diation (Hogan & Weynmann 1987; Gould & Weinberg 
1996) or on highly symmetric semi-analytical models for 
the gas distribution (Zheng & Miralda-Escude 2002b). 
Our results show that simple models (e.g. Gould & 
Weinberg 1996) tend to overpredict the Lya flux emitted 
from optically thick clouds. In fact, we identified two 
effects that reduce the fluorescent Lya flux (and mod- 
ify the spectral energy distribution) with respect to the 
widespread static and plane-parallel model. 

Velocity effect - The velocity field inside and around 
the shielding layers of a gas cloud influences the emerg- 
ing line profile. The symmetry of the double humped 
spectrum is lost and, in most cases, one of the two peaks 
is severely suppressed. On average, the SB of a cloud is 
reduced by 25% with respect to the static situation. 

Geometric effect - For anisotropic illumination and in 
the presence of a strong ionizing flux, the thickness of 
the shielding layer is comparable to the size of the gas 
cloud. In this case, the angular distribution of the emerg- 
ing radiation is very different than in the plane-parallel 
approximation. For instance, close to a quasar, a cloud 
emits much less than predicted by the slab model in the 
direction of the quasar and much more in the other di- 



rections. 

The importance of these effects (in particular of the 
angular redistribution of Lya photons) depends on the 
intensity of the impinging radiation. In equation (|16|) . we 
provided a fitting function for the maximum Lya bright- 
ness of optically thick sources as a function of the incident 
ionizing flux. In Figure ITT1 we presented our predictions 
for the number density of fluorescent sources with differ- 
ent sizes. 

These results are consistent with the recent null detec- 
tion by Francis & Bland-Hawthorn (2004) and represent 
a fundamental tool for planning future observations. 

We would like to thank Piero Madau for suggest- 
ing to consider the problem of anisotropic illumina- 
tion in the presence of a quasar. We benefited from 
extensive discussions with Francesco Haardt, Martin 
Haehnelt, Tom Theuns and all the participants to the 
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cently held at the Kavli Institute for Theoretical Physics 
in Santa Barbara. We thank the anonymous referee for 
useful suggestions which improved the presentation of 
our work. CP and FM acknowledge support from the 
Zwicky Prize Fellowship program at ETH-Ziirich and SC 
from the Swiss National Science Foundation. 



APPENDIX 

RADIATIVE TRANSFER FOR A DIFFUSE BACKGROUND 

The fluorescent Lya flux from a source is proportional to the number of hydrogen ionizations taking place in the 
corresponding gas cloud. In order to simplify the radiative transfer problem for continuum radiation, in equation J2J 
we only followed the propagation of ionizing photons along the 6 principal directions of the simulation box. In this 
appendix, we first provide justification for the multiplicative factor 27r/3 which appears in equation @ and then test 
the 6-direction approximation. 

The factor 27r/3 in equation (j2J is obtained by discretizing the left-hand side in equation JQ) along 6 preferred 
directions. Note, however, than equation applies to an infinitesimal volume element while our discretized version is 
used to describe finite cubic cells. It is clear that the approximation is correct for optically thin cells, but where does 
it break down? Let us consider a homogeneous cube of side L (corresponding to an optical depth At co h) illuminated 
by a uniform and monochromatic background with intensity / (photons per unit time, surface and solid angle). The 
ionization rate within the cube generated by the radiation background impinging on one face can be written as 
a it / L 2 e~ ATcc11 , with a a numerical coefficient of order unity. We use a Monte Carlo method to compute a as a 
function of Ar cc n. Our results are plotted in Figure ITS! When the cell is optically thin a ~ 2/3, while a — > 1 when 
Ar ce n — ► oo. A rather sharp transition between the two asymptotic regimes is observed for Ar cc u ~ 10. Note that, in 
the presence of a uniform background, the optically thin approximation adopted in the main text is accurate even for 
Ar cc n = 1, where a = 0.70. Given that, in our adaptively refined cosmological simulation, ~ 93% of recombinations 
take place in cells with Ar cc u < 1, we used a = 2/3 in our calculations. 

All this, however, applies to an isolated cell exposed to the ionizing background. We now want to test how accurately 
equation J5J works in the shielding layer of an optically thick cloud. We thus consider a spherically symmetric cloud 
which is optically thin at his external boundary and thick at the center (logAr co n = 2 — 4(r/i£), with r the radial 
coordinate which vanishes at the cloud center and R the characteristic size of the cloud, which roughly mimics one of 
the clouds in our refined simulation). We then illuminate the outer boundaries of the cloud with a uniform (over the 
external 2ir sr) and monochromatic background. This is injected from the surface of a cube of side 2R centered on the 
cloud and divided into 151 3 mesh points. In Figure H*2l we compare the outcome of a detailed radiative transfer code 
(Porciani & Madau, in preparation) with our approximated method. Note that in our six-directions approximation 
the ionizing flux penetrates slightly deeper into the cloud with respect to the exact solution. This happens because 
our method ignores photon trajectories which are oblique to the principal axes of the box. On the other hand, we 
find that our approximation produces 92% of the ionizations that take place in the cloud. Given the simplicity of our 
algorithm, this is a remarkable achievement. Note that, adopting the optically thick approximation (a — 1), we would 
have overestimated the number of recombinations by 38%. These figures are rather stable and, for realistic cases, do 
not depend on the details of the optical depth distribution. 
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Fig. 12. — Left: The ionization rate generated by a uniform UV background (of intensity /) penetrating in a cubic cell (of size L) 
through one of its faces as a function of the optical depth variation within the cell. The solid line represents the ionization rate (in units 
of 7r 1 1? e _ATco11 ) while the dashed line marks the optically thin solution a = 2/3. Right: Radial profile of the ionization rate (arbitrary 
units) for the mock cloud described in Appendix [X] The exact solution (solid) is compared with our six-direction approximation in the 
optically thin (dashed) and optically thick (dotted) cases. 
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